/*
 *  RKModMaxwell.cpp
 *  ModMaxwell
 *
 *  Created by Michael Barriault on 10-10-01.
 *  Copyright 2010 MikBarr Studios. All rights reserved.
 *
 */

#include "RKModMaxwell.h"
#include <T3nsors/Others.h>
#include <T3nsors/Macros.h>
#include <cmath>
#include <cstdio>
using namespace T3;

#define Phi x.S(0)
#define Psi x.S(1)
#define E x.V(0)
#define B x.V(1)

#define dPhi2 dx2.S(0)
#define dPsi2 dx2.S(1)
#define dE2 dx2.V(0)
#define dB2 dx2.V(1)
#define dPhi dx.S(0)
#define dPsi dx.S(1)
#define dE dx.V(0)
#define dB dx.V(1)

const real zero = 1e-6;

Set RKModMaxwell::Dots(Set x) {
    Domain O = x.O;
    Set dx(O,x.N,x.s,x.v);
    Scalar divE = div(E);
    Vector EtB = E%B;
    Scalar EcE = E*curl(E);
    Scalar BcB = B*curl(B);
    //    Vector J = div(E)*(E%B)+1./(4*M_PI)*(E*curl(E)-B*curl(B))*B;
    Vector J(O,3,"J");
    PFORO FORa J(a)[o] = divE[o]*EtB(a)[o]+(BcB[o]-EcE[o])*B(a)[o];
    Scalar B2 = B*B;
    real B0 = mu/pow(R,3.);
    real B02 = pow(zero*B0,2);
    PFORO FORa if ( B2[o] >  B02 ) J(a)[o] /= 4*M_PI*B2[o]; else J(a)[o] = 0;
    Scalar EB = E*B;
    Vector EBEB = Vector(O,3);
    PFORO FORa EBEB(a)[o] = EB[o]*(E(a)[o]+B(a)[o]);
    dPhi = LinCom(-K, Phi, /*1., div(E),*/ -epsilon, diss(Phi)); // dPhi
    dPsi = LinCom(-K, Psi, 1., div(B), -epsilon, diss(Psi)); // dPsi
    dE = LinCom(1., curl(B), step(K), diff(Phi), -4*M_PI, J, -I, EBEB, -epsilon, diss(E)); // dE
    dB = LinCom(-1., curl(E), step(K), diff(Psi), -I, EBEB, -epsilon, diss(B)); // dB
    Boundaries(dx);
    return dx;
}

void RKModMaxwell::Mod(Set& x) {
    Domain O = x.O;
    // Boundary conditions
	real sint = pole;
	real cost = 1-pole;
    PFOR(j,O(1)) FOR(i,O(0)) {
        B(0)(i,j,0) = 0;
        B(1)(i,j,0) = mu*sint/(R*R*R);
        B(2)(i,j,0) = 2*mu*cost/(R*R*R);
        E(0)(i,j,0) = 0;
        E(1)(i,j,0) = -2*Omega*mu*sint*cost/(R*R);
        E(2)(i,j,0) = Omega*mu*sint*sint/(R*R);
    }
    Scalar E2 = E*E;
    Scalar B2 = B*B;
    Scalar EB = E*B;
    real E02 = pow(zero*Omega*mu/(R*R),2);
    real B02 = pow(zero*mu/(R*R*R),2);
    // Remove proj_B E
/*    Vector F(O,3);
    PFORO if ( B2[o] > B02 ) FORa F(a)[o] = EB[o]/B2[o]*B(a)[o];
    PFORO FORa E(a)[o] -= F(a)[o];*/
    // Shrink so E<B
    real f = 1;
//    FORO if ( E2[o] >= B2[o] ) { f = 0.9; break; }
//    if ( f<1 ) PFORO FORa E(a)[o] *= f;
}

const real RT2 = sqrt(2);
const real RT3 = sqrt(3);

#define PT(obj) obj(i,j,k)
#define SUM(obj,one) (obj(i,j,k) one(i,j,k))
#define NORM2(obj,one,two) (obj(i,j,k) + (one(i,j,k) two(i,j,k))/RT2)
#define NORM3(obj,one,two,thr) (obj(i,j,k) + (one(i,j,k) two(i,j,k) thr(i,j,k))/RT3)

void BoundariesXYZ(Set dx, Set& dx2, int i, int j, int k) {
    Domain O = dx.O;
    real pmx = PM(i,O(0));
    real pmy = PM(j,O(1));
    real pmz = PM(k,O(2));
    PT(dPhi2) = 0.5*NORM3(dPhi, +pmx*dE(0), +pmy*dE(1), +pmz*dE(2));
    PT(dPsi2) = 0.5*NORM3(dPsi, +pmx*dB(0), +pmy*dB(1), +pmz*dB(2));
    PT(dE2(0)) = 0.5*NORM3(dE(0), +pmx*dPhi, +pmy*dB(2), -pmz*dB(1));
    PT(dE2(1)) = 0.5*NORM3(dE(1), -pmx*dB(2), +pmy*dPhi, +pmz*dB(0));
    PT(dE2(2)) = 0.5*NORM3(dE(2), +pmx*dB(1), -pmy*dB(0), +pmz*dPhi);
    PT(dB2(0)) = 0.5*NORM3(dB(0), +pmx*dPsi, -pmy*dE(2), +pmz*dE(1));
    PT(dB2(1)) = 0.5*NORM3(dB(1), +pmx*dE(2), +pmy*dPsi, -pmz*dE(0));
    PT(dB2(2)) = 0.5*NORM3(dB(2), -pmx*dE(1), +pmy*dE(0), +pmz*dPsi);
}

void BoundariesXY(Set dx, Set& dx2, int i, int j, int k) {
    Domain O = dx.O;
    real pmx = PM(i,O(0));
    real pmy = PM(j,O(1));
    PT(dPhi2) = 0.5*NORM2(dPhi, +pmx*dE(0), +pmy*dE(1));
    PT(dPsi2) = 0.5*NORM2(dPsi, +pmx*dB(0), +pmy*dB(1));
    PT(dE2(0)) = 0.5*NORM2(dE(0), +pmx*dPhi, +pmy*dB(2));
    PT(dE2(1)) = 0.5*NORM2(dE(1), -pmx*dB(2), +pmy*dPhi);
    PT(dE2(2)) = 0.5*NORM2(dE(2), +pmx*dB(1), -pmy*dB(0));
    PT(dB2(0)) = 0.5*NORM2(dB(0), +pmx*dPsi, -pmy*dE(2));
    PT(dB2(1)) = 0.5*NORM2(dB(1), +pmx*dE(2), +pmy*dPsi);
    PT(dB2(2)) = 0.5*NORM2(dB(2), -pmx*dE(1), +pmy*dE(0));
}

void BoundariesYZ(Set dx, Set& dx2, int i, int j, int k) {
    Domain O = dx.O;
    real pmy = PM(j,O(1));
    real pmz = PM(k,O(2));
    PT(dPhi2) = 0.5*NORM2(dPhi, +pmy*dE(1), +pmz*dE(2));
    PT(dPsi2) = 0.5*NORM2(dPsi, +pmy*dB(1), +pmz*dB(2));
    PT(dE2(0)) = 0.5*NORM2(dE(0), +pmy*dB(2), -pmz*dB(1));
    PT(dE2(1)) = 0.5*NORM2(dE(1), +pmy*dPhi, +pmz*dB(0));
    PT(dE2(2)) = 0.5*NORM2(dE(2), -pmy*dB(0), +pmz*dPhi);
    PT(dB2(0)) = 0.5*NORM2(dB(0), -pmy*dE(2), +pmz*dE(1));
    PT(dB2(1)) = 0.5*NORM2(dB(1), +pmy*dPsi, -pmz*dE(0));
    PT(dB2(2)) = 0.5*NORM2(dB(2), +pmy*dE(0), +pmz*dPsi);
}

void BoundariesXZ(Set dx, Set& dx2, int i, int j, int k) {
    Domain O = dx.O;
    real pmx = PM(i,O(0));
    real pmz = PM(k,O(2));
    PT(dPhi2) = 0.5*NORM2(dPhi, +pmx*dE(0), +pmz*dE(2));
    PT(dPsi2) = 0.5*NORM2(dPsi, +pmx*dB(0), +pmz*dB(2));
    PT(dE2(0)) = 0.5*NORM2(dE(0), +pmx*dPhi, -pmz*dB(1));
    PT(dE2(1)) = 0.5*NORM2(dE(1), -pmx*dB(2), +pmz*dB(0));
    PT(dE2(2)) = 0.5*NORM2(dE(2), +pmx*dB(1), +pmz*dPhi);
    PT(dB2(0)) = 0.5*NORM2(dB(0), +pmx*dPsi, +pmz*dE(1));
    PT(dB2(1)) = 0.5*NORM2(dB(1), +pmx*dE(2), -pmz*dE(0));
    PT(dB2(2)) = 0.5*NORM2(dB(2), -pmx*dE(1), +pmz*dPsi);
}

void BoundariesX(Set dx, Set& dx2, int i, int j, int k) {
    Domain O = dx.O;
    real pm = PM(i,O(0));
    PT(dPhi2) = 0.5*SUM(dPhi,+pm*dE(0));
    PT(dPsi2) = 0.5*SUM(dPsi,+pm*dB(0));
    PT(dE2(0)) = 0.5*SUM(dE(0),+pm*dPhi);
    PT(dE2(1)) = 0.5*SUM(dE(1),-pm*dB(2));
    PT(dE2(2)) = 0.5*SUM(dE(2),+pm*dB(1));
    PT(dB2(0)) = 0.5*SUM(dB(0),+pm*dPsi);
    PT(dB2(1)) = 0.5*SUM(dB(1),+pm*dE(2));
    PT(dB2(2)) = 0.5*SUM(dB(2),-pm*dE(1));
}

void BoundariesY(Set dx, Set& dx2, int i, int j, int k) {
    Domain O = dx.O;
    real pm = PM(j,O(1));
    PT(dPhi2) = 0.5*SUM(dPhi,+pm*dE(1));
    PT(dPsi2) = 0.5*SUM(dPsi,+pm*dB(1));
    PT(dE2(0)) = 0.5*SUM(dE(0),+pm*dB(2));
    PT(dE2(1)) = 0.5*SUM(dE(1),+pm*dPhi);
    PT(dE2(2)) = 0.5*SUM(dE(2),-pm*dB(0));
    PT(dB2(0)) = 0.5*SUM(dB(0),-pm*dE(2));
    PT(dB2(1)) = 0.5*SUM(dB(1),+pm*dPsi);
    PT(dB2(2)) = 0.5*SUM(dB(2),+pm*dE(0));
}

void BoundariesZ(Set dx, Set& dx2, int i, int j, int k) {
    Domain O = dx.O;
    real pm = PM(k,O(2));
    PT(dPhi2) = 0.5*SUM(dPhi,+pm*dE(2));
    PT(dPsi2) = 0.5*SUM(dPsi,+pm*dB(2));
    PT(dE2(0)) = 0.5*SUM(dE(0),-pm*dB(1));
    PT(dE2(1)) = 0.5*SUM(dE(1),+pm*dB(0));
    PT(dE2(2)) = 0.5*SUM(dE(2),+pm*dPhi);
    PT(dB2(0)) = 0.5*SUM(dB(0),+pm*dE(1));
    PT(dB2(1)) = 0.5*SUM(dB(1),-pm*dE(0));
    PT(dB2(2)) = 0.5*SUM(dB(2),+pm*dPsi);
}

void BoundariesXYPeriodic(Set dx, Set& dx2) {
	Domain O = dx.O;
	PFOR(k,O(2)) {
		FOR(i,O(0)) {
			dPhi2(i,0,k) = dPhi2(i,O(1)-1,k) = 0.5*( dPhi(i,0,k) + dPhi(i,O(1)-1,k) );
			dPsi2(i,0,k) = dPsi2(i,O(1)-1,k) = 0.5*( dPsi(i,0,k) + dPhi(i,O(1)-1,k) );
			FORa dE2(a)(i,0,k) = dE2(a)(i,O(1)-1,k) = 0.5*( dE(a)(i,0,k) + dE(a)(i,O(1)-1,k) );
			FORa dB2(a)(i,0,k) = dB2(a)(i,O(1)-1,k) = 0.5*( dB(a)(i,0,k) + dB(a)(i,O(1)-1,k) );
		}
		FOR(j,O(1)) {
			dPhi2(0,j,k) = dPhi2(O(0)-1,j,k) = 0.5*( dPhi(0,j,k) + dPhi(O(0)-1,j,k) );
			dPsi2(0,j,k) = dPsi2(O(0)-1,j,k) = 0.5*( dPsi(0,j,k) + dPsi(O(0)-1,j,k) );
			FORa dE2(a)(0,j,k) = dE2(a)(O(0)-1,j,k) = 0.5*( dE(a)(0,j,k) + dE(a)(O(0)-1,j,k) );
			FORa dB2(a)(0,j,k) = dB2(a)(O(0)-1,j,k) = 0.5*( dB(a)(0,j,k) + dB(a)(O(0)-1,j,k) );
		}
		dPhi2(0,0,k) = dPhi2(0,O(1)-1,k) = dPhi2(O(0)-1,0,k) = dPhi2(O(0)-1,O(1)-1,k) = 0.25*( dPhi(0,0,k) + dPhi(0,O(1)-1,k) + dPhi(O(0)-1,0,k) + dPhi(O(0)-1,O(1)-1,k) );
		dPsi2(0,0,k) = dPsi2(0,O(1)-1,k) = dPsi2(O(0)-1,0,k) = dPsi2(O(0)-1,O(1)-1,k) = 0.25*( dPsi(0,0,k) + dPsi(0,O(1)-1,k) + dPsi(O(0)-1,0,k) + dPsi(O(0)-1,O(1)-1,k) );
		FORa dE2(a)(0,0,k) = dE2(a)(0,O(1)-1,k) = dE2(a)(O(0)-1,0,k) = dE2(a)(O(0)-1,O(1)-1,k) = 0.25*( dE(a)(0,0,k) + dE(a)(0,O(1)-1,k) + dE(a)(O(0)-1,0,k) + dE(a)(O(0)-1,O(1)-1,k) );
		FORa dB2(a)(0,0,k) = dB2(a)(0,O(1)-1,k) = dB2(a)(O(0)-1,0,k) = dB2(a)(O(0)-1,O(1)-1,k) = 0.25*( dB(a)(0,0,k) + dB(a)(0,O(1)-1,k) + dB(a)(O(0)-1,0,k) + dB(a)(O(0)-1,O(1)-1,k) );
	}
}

void RKModMaxwell::Boundaries(Set& dx) {
    Set dx2 = dx;
    Domain O = dx.O;
    BoundariesXYPeriodic(dx, dx2);
    PFORO {
        int i = o % O(0);
        int j = ((o - i)/O(0)) % O(1);
        int k = ((o - i - j*O(0))/(O(0)*O(1))) % O(2);
        
        /*if ( (i == 0 or i == O(0)-1 ) and ( j == 0 or j == O(1)-1 ) and ( k == 0 or k == O(2)-1 ) )
            BoundariesXYZ(dx, dx2, i, j, k);
        else if ( (i == 0 or i == O(0)-1 ) and ( j == 0 or j == O(1)-1 ) )
            BoundariesXY(dx, dx2, i, j, k);
        else if ( (i == 0 or i == O(0)-1 ) and ( k == 0 or k == O(2)-1 ) )
            BoundariesXZ(dx, dx2, i, j, k);
        else if ( (j == 0 or j == O(1)-1 ) and ( k == 0 or k == O(2)-1 ) )
            BoundariesYZ(dx, dx2, i, j, k);
        else if ( i == 0 or i == O(0)-1 )
            BoundariesX(dx, dx2, i, j, k);
        else if ( j == 0 or j == O(1)-1 )
            BoundariesY(dx, dx2, i, j, k);
        else*/ if ( k == 0 or k == O(2)-1 )
            BoundariesZ(dx, dx2, i, j, k);
    }
    dx = dx2;
}



